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I.  INTRODUCTION 


In  recent  decades  z-transform  techniques  have  found  use  in  the  digital 
simulation  of  continuous  linear  systems.  One  of  these  techniques^  is  compared 
with  a  classical  method^  and  possible  methods  for  improvement  of  the  technique 
are  suggested.  The  emphasis  will  be  in  making  the  time  step,  T,  as  large  as 
possible  for  speed  in  real-time  simulation  or  economy  in  Monte  Carlo  studies. 
Another  possibility  would  be  improvements  in  the  fidelity  of  the  simulation. 
The  approach  would  also  find  application  in  microprocessor  and  microcomputer 
programming . 

Since  the  steady  state  response  of  a  system  is  of  interest,  the  response 
of  the. plant  to  the  Heaviside  unit  step, 

10,  t  <  0 

,  (1) 

l.tiO 

will  be  analyzed.  Other  inputs  considered  are  the  damped  exponential  and 
sine-wave . 

The  single  pole  filter  will  serve  as  a  test  case  since  it  is  the  basic 
transfer  function.  The  technique  implies  the  use  of  the  Jordan  canonical 
form.^ 

Before  proceeding  further,  some  commei 
on  the  conventions  used  in  this  report  are 
warranted  (Figure  1).  From  the  "sifting" 
property  of  the  Dirac  delta  distribution 

oo 

f (nT)  =  f  f(t)  6(nT  -  t)  dt  . 

— OO 

The  Laplace  transform  is  defined  to  be 
00 

F(s)  =  T  f(t)  u(t)  e  st  dt 

_on 


f(t)  f(nT) 

- c r  T  o - 

F(s)  F(z) 

Figure  1.  Time,  Sample  Data, 
Frequency  and  Sequency 
"Domains. " 

(2) 


(3) 


and  the  z-transform,  a  "discrete"  Laplace  transform. 


F(z)  -  f  Z  f<0  «(nT  -  t)  u(t)  e'St  dt  ,  (4) 

—CO  n  =  —CD 

which  readily  reduces  to 


F(z) 


Z  f (nT)  u(nT)  e  snT  . 

n  =  -oo 


(5) 
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PRECEDING  PACE  BLANK-NOT  flLt«£D 


For 


-sT 

z  =  e  , 

Equation  (5)  becomes^ 


F(z)  =  £  f(nT)  zn  . 
n=0 

The  modified  z-transform  is 


(6) 


F(z,a)  =  I  F(nT  +  aT)  zn  ,  (7) 

n=0 


and  finally 


L,[f (n)  (t)l  =  snF(s)  -  t  s11"4"1  f&)(0)  ,  (8) 

Jl=0 

which  is  used  to  incorporate  non-zero  initial  conditions. 

Since  the  z-transform  is  developed  from  the  known  time  domain  solution  of 
a  linear  differential  equation,  the  recurrences  developed  for  a  given  input 
and  plant  are  exact  for  any  size  time  step,  T.'4 

Transforming  the  linear  differential  equation  for  the  single  pole  filter 
with  a  unit  step  input , 

y(t)  +  a  y (t)  =  u(t)  ,  (9) 

to  the  Laplace  domain,  using  Equation  (8)  yields 


[s  Y(s)  +  y(0)]  +  a  Y(s)  =  . 

Collecting  terms  yields 


(10) 


Y(s) 


s  +  a  s(s  +  a) 

Taking  the  z-transform  of  Equation  (11), 


(11) 


y«  -  y(o)  zfpLj]  +  '  <12) 

Since  the  initial  condition  is  a  constant  in  the  Laplace  domain  (an  impulse  in 
the  time  domain),  the  Raggazzini-Zadeh  identity  applies. 


4 


*  ✓ 


(20c) 


y[n  -  (1  -  a)]  =  1  -  a[y(n  -  1)  +  aT  y(n  -  1)]  . 

Substituting  and  collecting  terms, 

y(n)  a  (1  -  aT  +  (aT)2/2)  y(n  -  1) 

+  [l  -  (1  -  aT  +  (aT)2/2)]/a  ;  (21) 

the  a's  cancel  out! 

Note  that 


e  aT  -  1  -  aT  +  (aT)2/2  (22) 

is  the  (2/0)  Pade  approximation  for  the  exponential.^  This  is  not  unexpected 
since  the  Runge-Kutta  integrators  were  developed  by  truncating  the  Taylor 
series,  and  the  (N/0)  Pade  approximations  for  the  exponential  are  just  the 
truncated  Taylor  series. 

In  the  sequency  domain.  Equation  (14)  would  be 

Y(z)  =  V(0)  +  fl  -  (1  -  aT  +  (aT)2/2)~]z/a(l  -  z)  (23) 

1  -  (1  -  aT  +  (aT)2/2)z 

To  develop  a  recurrence  for  the  error7  of  the  Runge-Kutta  integrator, 
subtract  Equation  (4)  from  Equation  (23): 

e(z)  =  y(0)  +  fl  -  (1  -  aT  +  (aT)2/2)1 z/a(l  -  z) 

1  -  (1  -  aT  +  (aT)2/2)z 


-aT. 


_  y(0)  +  (1  -  e  )z/a(l  -  z) 


-aT 
e  z 


(24) 


Without  proceeding  further,  it  is  readily  apparent  that  the  initial  condition 
does  not  cancel  out. 


In  Figure  2,  the  (2/0)  Pade 
approximation  is  plotted  with  the 
exact  function.  The  approximation 
not  only  becomes  inaccurate  when  aT 
is  too  large,  but  exceeds  1  for  aT 
greater  than  2. 


Figure  2. 


(2/0)  Pade  for  e 


-aT 
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III.  A  Z -TRANSFORM  APPROACH 


Of  course,  che  z-transform  is  exact  when  the  input  and  plant  are  both 
given  a  priori.  In  a  digital  simulation,  the  input  is  a  sequence  of  values 
and  not  a  known  function,  that  is,  it  is  not  generally  known  a  priori. 


In  this  case,  approximations  for  the  z-transform  of  a  product  in  the 
Laplace  domain  are  required.  The  problem  is  to  approximate  the  convolution 
integral ^  in 


z[G(s)  F ( s )]  =  £  zn  u(nT)  f  g(t>  u(t)  f(nT  -  t)  u(nT  -  t)dt 

(25a) 

00  OO  kT+T 

=  l  Zn  u(nT)  I  /  g(t)  u(t)  f (nT  -  t)  u(nT  -  t)dt 

n  =  -oo  k  =  -°°  kT 

(25b) 

15  8 

Those  approximations  previously  developed  ’  *  are  shown  in  Table  1. 


The  second  order  Runge-Kutta  integrator  family  [see  Equation  (20a)]  leads  to 

kf  T  g(t)  e(nT  -  t)dt  *  T^l  -  g(kT)  f(nT  -  kT) 
kT 

+  g(kT  +  T)  f (nT  -  (k  +  a)T)  .  (26) 


Substituting  Equation  (26)  into  Equation  (25b)  would  lead  to  the  desired 
results  after  some  manipulation  of  the  series. 5  A  shortcut  is  available; 
note  that  on  the  right-hand  side  of  Equation  (26)  the  first  term  leads  to  Eular 
Convolution  and  the  second  to  Mean  Value  Convolution.  With  this  insight,  t  e 
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results  may  be  written  down  by  inspection  as  was  done  before  for  some  higher 
order  Runge-Kutta  convolutions , 5 

z[G(s)  F(s)J  -  T^l  -  G(z)  [f(z)  -  f(0)] 

+  G(z,a)[F(z,-a)  -  f(-aT)]  ,  (27) 

R-K(2,ot)  Convolution.  Eular  Convolution  and  Mean  Value  Convolution  are  just 
special  cases  of  R-K(2,a)  Convolution  (see  Table  2). 


TABLE  2.  SPECIAL  CASES  OF  R-K(2,ct)C 


For  a  single  real  pole  filter. 


FW  ■  7TT  ■ 

and  any  input,  G(s),  the  approximation  using  R-K(2,a)C  is 

*  i1  -  £)  -  *] 

+  ih)  ■  “T] 


r(F  -  h)  e'aTz  GW  +  T( 


-(1  -o)aT 


z  G(z,a) 


,  -aT 
1  -  e  z 


(29a) 


(29b) 


The  initial  condition  would  be  incorporated  exactly  as  in  Equations  (9)  -(16a). 
Equating  coefficients  of  like  powers  of  "z,"  one  has  the  recurrence 


y(n)  -  e  aT  y(n  -  1)  +  T^l  -  e  aT  g(n  -  1) 
+  (^i)e  ^  "°^aT  gCn  -  (1  -  o))j  ,  n  >  0  . 
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For  a  Heaviside  unit  step  input,  G(s)  =  —  » 


J(,  l\ -aT  /l\  -(l-a)aTl 

,[(i\(  i  \]  J 

LVsAs  +  a/J"  (1  -  Z)(l  -  e-aTz) 


While  the  denominator  of  Equation  (31)  is  exact,  in  the  numerator 


e"aT  =  1  -  aT 


K1  -  h) 


-aT  u. 
e  + 


-(1  -  a)aT' 


As  aT  becomes  very  large,  the  right-hand  side  of  Equation  (32)  approaches  unity. 
Plots  of  Equation  (32)  are  shown  in  Figure  3. 


Figure  3.  Plots  of  Equation  (32). 


A  single  integrator  will  integrate  a  Heaviside  unit  step  exactly  for  any 
size  time  step,^  but  the  recurrence  for  a  single  real  pole  filter  will  not. 
For  what  input,  if  any,  is  the  recurrence  exact? 

Two  exact  z-transforms  of  interest  are"* 


f _ 1 _ ]  =  _ (e~aT  -  e~bT)z _ 

l(s  +  b)(s  +  a)J  "  (b  _  a)(1  _  e-bTz)(i  _  e-aTz) 


(33a) 


T  e_a^z 


/ 1  — sT  %  * 

(1  -  e  z) 


(33b) 


*  / 


For  R-K(2,a)C,  Equation  (27),  one  has 


{(.  +  X  l  a)]  •  4 ' 2a)  (j  _  x  _  e-«Tz  -  1 


aaT 


■t  ~aT 

1  -  e  z 


aaT 

a 


(34a) 


T  z 


f(x  -  j)  ■~aT  +  (j) 


-aT  ,  /  1  \  -abT  -(1  -a)aT 
e  +  ( -rr  /  e  e 


-bT  w,  -aT  . 
(1  -  e  z)  (1  -  e  z) 


.  (34b) 


No  choice  of  a  leads  to  Equation  (33a).  For  a  =  b,  Equation  (34b)  becomes 


T  ~aT, 

X  e  z 


(1  -  e-aTz) (1  -  e~aTz) 


(33b) 


which  is  exact  for  any  aT!  This  implies  that  the  recurrence  for  a  single  real 
pole  filter  can  filter  a  decaying  exponential  exactly  for  any  size  time  step, 

T,  where  they  have  identical  time  constants  (a  =  b) . 

The  unit  step  into  a  single  integrator  is  just  a  special  case  (a  =  b  =  0)  . 
From  this  point  of  view,  the  single  integrator  is  just  a  single  pole  filter  at 
the  origin;  the  recurrence  for  the  single  pole  filter  reduces  to  the  classical 
integrator (s)  for  a  =  0. 


IV.  ADJUSTMENT  OF  UNIT  STEP  RESPONSE  VIA  RESIDUES 

It  is  generally  more  desirable  to  have  the  recurrence  for  a  single  real 
pole  filter  respond  exactly  to  the  Heaviside  unit  step.  To  determine  how  this 
might  be  accomplished,  it  is  necessary  to  first  determine  the  error  for  a  unit 
step  input  into  a  single  real  pole  filter.  To  develop  a  recurrence  for  the 
error, '  subtract  the  exact  z-transform  from  the  convolution  approximation; 
Equation  (31)  minus  Equation  (13), 

'(.>  -  C"  -  a  ,  os) 

(l-z)(l-eaz) 


where 


(36) 


Since  the  initial  condition  would 


be  incorporated  exactly. 


1,5 


it  would  vanish. 
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Equating  coefficients  of  like  powers  of  z, 
e(0)  =  0  , 

e(nT)  =  e~aT  e(nT  -  T)  +  [n  -  (1  -  e-aT)/a],  n  >  0  . 
Writing  a  few  terms 

e(T)  =  N  -  (1  -  e~aT) / a 
e(2T)  =  (1  +  e"aT)[N  -  (1  -  e“aT)/a]  , 
it  becomes  apparent  that 

e(nT)  =  (""I  e_kaT\  [fj  _  (l  -  e_aT)  /a]  . 

\k=0  / 


Factoring  yields 

1  -  e 


e(nT)  = 


/"f1  -kaT\ 

N  -  1 

Uoe  ) 

.(1  -  e-aT) / a 

To  remove  the  error,  introduce  a  residue  into  Equation  (29a), 

z[G<s)(rfi)]  ■ Res  T[(1  -  n)  g<2)  e“T 

+  {ja')  6_(1  _a)aTz  G(z,o)J/(l  -  e“aTz)  . 

Equation  (39)  would  then  become 


(ny-  -kaT\ 

Res  N  1 

Uoe  ) 

.(1  -  e"aT)/a  . 

Setting  the  quantity  in  the  brackets  equal  zero  gives 


Res 


(1  -  e~aT) / a 


For  this  residue,  Equation  (42),  the  recurrence 


(37a) 

(37b) 

(38a) 

(38c) 

(39) 


(40) 

(41) 

(42) 
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y(n) 


-  e  aT  y (n  -  1)  +  Res  T ^1  -  e  aT  g(n  -  1) 
+  e  ^  ~  a^aT  g(n  -  (1  -  a))J  ,  n  >  0 


(43) 


would  be  exact  for  a  unit  step  input  for  any  size  time  step,  T! 

O 

For  a  =  1,  that  is,  Trapezoidal  Convolutions 

N  =  ~  (1  +  e_aT)  (44) 

and 


Res  = 


(1  -  e~aT) / a 
i  (1  +  e~aT)  ’ 


Res  =  ^  tanh  (^)  . 


(45a) 


(45b) 


In  this  case,  the  recurrence.  Equation  (43)  ,  would  be 

tanh 


y(n)  a  e  aT  y(n  -  1)  + 


[e-a  g(n  -  1)  +  g(n)]  ,  n  >  0  , 

(46) 


instead  of 


y(n)  =  e  aT  y(n  -  1)  +  -y-  [e  31  g(n  -  1)  +  g(n)J  ,  n  >  0  , 


X  r  -aT 
2 


(47) 


V.  ACCURACY  ANALYSIS 

Instead  of  the  error  for  a  unit  step  input,  consider  the  error  for  an 
exponential  input  into  a  single  pole  filter.  Equations  (35b)  minus  Equation 
(33)  yields 


e(z) 


z  N  -  z(e~aT  -  e  ^T)/(b  -  a) 

-bT  . ,,  -aT  . 

(1  -  e  z)(l  -  e  z) 


(48) 


Since  the  initial  condition  is  incorporated  exactly  in  either  case,  it 
vanishes. 


12 


*  ✓ 


Dividing  by  the  exact  z-transform,  Equation  (33a), 


Z 


e(z) 


.(£ 


+  b)  (s  +  a) 


(e 


-aT 


_ N _ 

e  bT)/(b  -  a) 


-  1  . 


(49) 


The  rationale  for  Equation  (49)  is  as  follows:  Since  the  initial  condi¬ 
tion  is  incorporated  exactly,  any  error  due  to  the  approximations  is  in  the 
lower  channel  of  Figure  4.  Dividing  by  Y(z)  would  introduce  the  initial 
condition  into  the  expression  for  the  (percent)  error. 


Figure  4.  Single  Real  Pole  Filter. 


Only  the  error  in  the  lower  channel  is  being  considered  and  not  the  error 
in  Y(z).  Of  course,  if  the  initial  conditions  are  zero,  there  is  no  diffi¬ 
culty  because  there  is  no  difference.  In  this  case, 

Y(z)  *  Z[(S  +  M(,  +  «)]  •  <50) 

There  is  ample  precedence  for  this  approach  in  transfer  function  analysis. 
To  have  the  ratio  of  output  to  input  all  the  initial  conditions  must  be  zero. 

aT 

Factoring  out  e  , 


(e 


-aT 


N _ 

e"bT)/(b  -  a) 


(b  -  a)T 
2a 


-a(b  -  a)T 


+  (2q  -  1) 


-(b-a)T 
e  -  1 


(51) 


In  the  preceding  it  has  been  tacitly  implied  that  "a"  and  "b"  are  both 
real  but  there  is  nothing  to  prevent  them  from  being  imaginary  or  complex. 

In  fact,  this  approach  was  motivated  by  an  analysis  of  a  sine  wave  input 

(b  =  oj)  into  single  integrators  (a  =  0).3,9,10  Actually,  the  input  considered 

was 


IojT 

e 


cos  wT  +  i  sin  u>T  , 


which  yields  both  amplitude  and  phase  response. 
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(52) 
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For 


b  =  a  +  io) 


and 


a  =  E  +  ift 

the  right-hand  side  of  Equation  (51)  becomes 


[(a  -  E)  +  i(u»  -  fl)1  T  (e-at(q-^)  +  K->-«)lT  ±  (2q  -  1)) 

2a(e-t(a_5:)  +  1(w_fi>lT  -  l) 


(52) 

(53) 


(55) 


Substituting  Eular's  Rule,  Equation  (52),  multiplying  top  and  bottom  by  the 
complex  conjugate  of  the  denominator  and  collecting  term:  for  the  real  part 


Re 


(o  -  E)T  A  -  (ai  -  fl)T  B 
2  a  D 


(56) 


and  for  the  imaginary  part 


Im 


i(oi  -  fl)T  A  +  i(g  -  E)T  B 
2  a  D 


(57) 


where 

D  =  e-2(<J-Z)T  _  2e-(a-E)T  cos  _  fi)T]  +  1  (58) 

A  =  e-d  +  «)(cr-  E)T  cos[-(1  _  a)(w  _  n)T] 

+  (2a  -  1)  e"(0"S)T  cos[(w-R)T]  -  1 

-  e-a^a  cos[a(di  -  8)t3  (59) 

B  =  e~d  +  a)  (a  -  E)T  sinj-(1  _  a)  (u)  _n)T] 

+  (2a  -  1)  e"(0_J:>T  sin  [(a:  -  fi)T] 

+  e  a^°  sin[a(a>  -  0)t3  •  (60) 

Of  course,  to  compute  the  amplitude  ratio: 

(Re^  +  Im^)*5  (61) 


while  the  phase  error  is 

arc  tan  (Im/Re)  .  (62) 

Note  that  if  either  w  »  fl  or  o  ■  I,  there  is  no  phase  error  when  a  *  1. 
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The  ranges  of  values  in  Figures  5-18  require  some  explanation.  The 
imaginary  part,  (w  -  ft)T,  ranges  from  0  to  ±ir(±180°).  The  upper  value  is  the 
Shannon  sampling  limit  of  two  samples  per  cycle.  The  curves  for  ±ir/2 (±90°) , 
that  is,  twice  the  Shannon  limit  or  four  samples  per  cycle,  are  a  more 
realistic  boundary  on  performance.  The  real  part,  (a  -  l)T,  ranges  from  -4 
to  +4.  A  realistic  working  limit  appears  to  be  ±2.  If  the  realistic  limit 
had  been  used  on  the  plots,  it  would  not  have  been  clear  what  happens  when 
these  limits  are  exceeded. 


VI.  RESULTS  AND  RECOMMENDATIONS 

A  "convolution"  approximation  based  upon  the  Runge-Kutta  second  order 
integrators  was  derived.  Equation  (27).  For  a  unit  step  into  a  single  real 
pole  filter  the  approximation,  R-K(2,ci)C,  proved  superior  to  the  Runge-Kutta 
integrators,  R-K(2,a)I;  compare  Figure  2  with  Figure  3. 

It  was  then  shown  that  when  the  input  and  the  plant  are  identical  single 
poles,  the  R-K(2,a)C  approximation  is  exact  for  any  a.  Equation  (34).  This 
rather  interesting  observation  was  then  used  to  develop  a  residue,  Equation  (42), 
which  makes  the  approximation  exact  for  a  unit  step  input  to  a  single  real 
pole  filter. 

Then  an  analysis  was  made  of  the  amplitude  and  phase  response  of  a  complex 
exponential  into  a  complex  pole  for  R-K(2,a)C.  A  value  of  a  =  2/3  offers 
excellent  amplitude  (Figure  13)  and  phase  (Figure  14)  response  over  a  wide 
range  of  complex  values. 

Since  the  amplitude  (Figure  5)  and  phase  (Figure  6)  response  for  R-K(2,1)C, 
Trapezoidal  Convolution,  is  not  so  well  behaved,  R-K(2,2/3)C  should  be  consid¬ 
ered  for  real  time  applications  when  a  large  time  step,  T,  is  required.  In  a 
"classical"  simulation  the  Runge-Kutta  second  order  integration  with  a  =  2/3 
is  recommended . 

There  may  be  some  difficulties  when  the  problem  is  implicit.  In  this 
case  Eular  Convolution  (a  =  0)  is  required  to  "predict"  the  advanced  value  for 
the  recurrence.  Van  der  Pol's  equation  would  serve  as  a  test  bed  for  the 
ability  of  the  technique  to  handle  implicit,  nonlinear  problems,  and  such  a 
study  has  been  proposed  as  a  follow-on  to  this  effort. 

This  technique  requires  further  development  and  refinement,  and  is  in  no 
way  complete  at  this  time,  but  results  to  date  are  very  promising. 
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Figure  17.  Amplitude  Ratio 
for  R-K(2,0.5)C. 


Figure  18.  Phase  Error  for 
R-K(2,0.5)C. 
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